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Abstract 

The general conditions under which the quadratic, uniform and monotonic 
convergence in the quasihnearization method of solving nonlinear ordinary 
differential equations could be proved are formulated and elaborated. The 
generalization of the proof to partial differential equations is straight forward. 
The method, whose mathematical basis in physics was discussed recently by 
one of the present authors (VBM) , approximates the solution of a nonlinear 
differential equation by treating the nonlinear terms as a perturbation about 
the linear ones, and unlike perturbation theories is not based on the existence 
of some kind of a small parameter. 

It is shown that the quasilinearization method gives excellent results when 
applied to different nonlinear ordinary differential equations in physics, such 
as the Blasius, Duffing, Lane-Emden and Thomas-Fermi equations. The first 
few quasilinear iterations already provide extremely accurate and numerically 
stable answers. 
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I. INTRODUCTION 



In a series of recent papers, [^0] the possibility of applying a very powerful approxima- 
tion technique called the quasilinearization method (QLM) to physical problems has been 
discussed. The QLM is designed to confront the nonlinear aspects of physical processes. The 
method, whose iterations are constructed to yield rapid convergence and often monotonicity, 
was originally introduced forty years ago by Bellman and Kalaba as a generalization 
of the Newton-Raphson method |^,^ to solve individual or systems of nonlinear ordinary 
and partial differential equations. Modern developments and applications of the QLM to 
different fields are given in a monograph 0. 

However, the QLM was never systematically studied or extensively applied in physics, 
although references to it can be found in well known monographs dealing with the 
variable phase approach to potential scattering, as well as in a few scattered research papers 
|[T0Hl3[. The reason for the sparse use of the QLM in Physics is that the convergence of 
the method has been proven only under rather restrictive conditions [0,0, which generally 
are not fulfilled in physical applications. Recently, though, it was shown [|l| by one of the 
present authors (VBM) that a different proof of the convergence can be provided which we 
will generalize and elaborate here so that the applicability of the method is extended to 
incorporate realistic physical conditions of forces defined on infinite intervals with possible 
singularities at certain points. 

In the first paper of the series [0], the quasilinearization approach was applied to the 
nonlinear Calogero equation in a variable phase approach to quantum mechanics and the 
results were compared with those of perturbation theory and the exact solutions. It was 
found analytically and by examples that the n-th approximation of the QLM exactly sums 
2" — 1 terms of perturbation theory. In addition, a similar number of terms is reproduced 
approximately. The number of the exactly reproduced perturbation terms thus doubles with 
each subsequent QLM approximation, and reaches, for example, 127 terms in the 6-th QLM 
approximation, 8191 terms in the 12-th QLM approximation, and so on. 
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The computational approach in the work [||] was mostly analytical, and therefore one 
was able to compute only two to three QLM iterations, mainly for power potentials. Only in 
the case of the potential, could the calculation of QLM iterations be done analytically 
for any n. 

The goal of the next work was, by dropping the restriction of analytical computation, 
to calculate higher iterations as well as to extend the analysis to non-power potentials, in 
order to better assess the applicability of the method and of its numerical stability and the 
convergence pattern of the QLM iterations. It was shown that the first few iterations already 
provide very accurate and numerically stable answers for small and intermediate values of 
the coupling constant and that the number of iterations necessary to reach a given precision 
only moderately increases for larger values of the coupling. The method provided accurate 
and stable answers for any coupling strengths, including for super singular potentials for 
which each term of the perturbation theory diverges and the perturbation expansion does 
not exist even for a very small coupling. 

The quasilinearization approach is applicable to a general nonlinear ordinary or partial n - 
th order differential equation in iV-dimensional space. In this paper, we consider the case of 
nonlinear ordinary differential equations in one variable which, unlike the nonlinear Calogero 
equation [§] considered in references contain not only quadratic nonlinear terms but 

various other forms of nonlinearity and not only a first, but also higher derivatives. Namely, 
we apply it to a panopoly of well-known and difficult nonlinear ordinary first, second and 
third order differential equations and show that again with just a small number of iterations 
one can obtained fast convergent and uniformly excellent and stable numerical results. 

The paper is arranged as follows: in the second chapter we present the main features of 
the quasilinearization approach, while in the third chapter we consider, as a warm-up exer- 
cise, a simple first-order differential equation with a nonlinear n-t\i power term and compare 
its exact analytic solution with the perturbation theory and with the QLM iterations in 
order to demonstrate the main features of the quasilinearization approach. In the next four 
chapters, we apply our method to four well known nonlinear ordinary second and third or- 
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der differential equations, namely to the Lane-Emden, Thomas- Fermi, Duffing, and Blasius 
equations, respectively. These equations have been extensively studied in the literature. For 



The results, convergence patterns, numerical stability, advantages of the method and its 
possible future applications are discussed in the final chapter. 



The aim of the quasilinearization method (QLM) of Bellman and Kalaba based 
on the Newton-Raphson method is to solve a nonlinear n-th order ordinary or partial 
differential equation in dimensions as a limit of a sequence of linear differential equations. 
This goal is easily understandable since there is no useful technique for obtaining the general 
solution of a nonlinear equation in terms of a finite set of particular solutions, in contrast 
to a linear equation which can often be solved analytically or numerically in a convenient 
fashion using superposition. In addition, the QL sequence should be constructed to assure 
quadratic convergence and, if possible, monotonicity. 

As we have mentioned in the Introduction, we will follow here the derivation outlined in 
ref . , which is not based, unlike the derivations in refs. , on a smallness of the interval 
and on the boundness of the nonlinear term and its functional derivatives, the conditions 
which usually are not fulfilled in physical applications. 

For simplicity, we limit our discussion to nonlinear ordinary differential equation in one 
variable on the interval [0, b], which could be infinite: 



a sample of recent papers see refs. [p!5|-|2^ and references therein. 



II. THE QUASILINEARIZATION METHOD (QLM) 




(2.1) 



with n boundary conditions 




(2.2) 



and 
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(6)) =0, k = l + l,...,n. (2.3) 

Here L*^"^ is linear n-th order ordinary differential operator and / and gi,g2, ,gn are 

nonlinear functions of u{x) and its — 1 derivatives u^'^^x), s = 1, ...n — 1. The more general 
case of partial differential equations in N- dimensional space could be considered in exactly 
the same fashion by changing the definition of L^"-* to be a linear n-th order differential 
operator in partial derivatives and x to be an N-dimensional coordinate array. 

The QLM prescription determines the r + 1-th iterative approximation Ur+i{x) to 

the solution of Eq. ( p.l|) as a solution of the linear differential equation 



L("V+i(x) = fiUr{x),ul'\x), ,U^;'-'\X),X) 

+ E(45i(a:) - ui^'Kx)) U^iurix),ul}\x), (x),a;), (2.4) 



n-1 



where u^^\x) = Ur{x), with linearized two-point boundary conditions 

j:{ul%{0)-ui^\0))g,^,4ur{0),ui'\0), ,u^r'\0),0) = 0, k = 1,.J (2.5) 



s=Q 

and 

n-1 



J2{ulfl,{b)-ul^\b))g,^,s,{ur{b),u^;\b), ,ui'"'\b),b) = 0, k = l + l,...,n. (2.6) 

s=0 



Here the functions f^(s) = df/du^^^ and g^ui^) = dgk/du^^\ s = 0,1,. ..,?7, — 
1 are functional derivatives of the functionals f{u{x),u^^\x),....u^^~^''{x),x) and 
gk{u{x),u^^\x), ....u^''^~^\x), x), respectively. []. 

The zeroth approximation Uq{x) is chosen from mathematical or physical considerations. 



^For example, in case of a simple nonlinear boundary condition u'{h)u{h) = c where c is a con- 
stant, one has g{r) = g{u{r),u' {r),r) = u'{r)u{r) so that Qu = u'{r) and g^' = u{r). The lin- 
earized boundary condition ^]6| has a form (ur+i(6) — Ur{b))u'^{b) + {u'j.^i{b) — u'^{h))u{b) = 
or {urJ^i(h)ur{h))' = {ur{b)ur{b)y SO the nonlinear boundary condition for the initial guess 
uo(6)uq(6) = c will be propagated to the linear boundary condition for the next iterations. 
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To prove that the above procedure yields a quadratic and often monotonic convergence 



s=0 



to the solution of Eq. |2.1| with the boundary conditions |2.2| and |2.3| , we follow reference |1| 
and consider a differential equation for the difference 6ur+i{x) = Ur+i{x) — Ur{x) between 
two subsequent iterations: 

+ fuMMx),ui^\x), 

- 6ulf\x) f^(s)U(^r - l{x),u^li{x), ,u';!'Si^\x),x)]. (2.7) 

The boundary conditions are similarly given by the difference of Eqs. |2.5| and |2.tj| for two 
subsequent iterations: 

n-l 

Y,[Suli,{o) <7,„,.)K(o),w«(o), ,u^r'\o),o) 

-6ui'\0) (7fc„(.)K_i(0),M(i\(0), ,ut~i'\0),0)] = 0, 

k = l,...l (2.8) 



and 



E l^^i (b) {ur (b) , (6) , , ^l'^- ^) (&) , fe) 

s=0 

-5w(^)(6) <7,„(.)(«,_i(6),t.;i\(6), ,4ri^)(^),&)] = 0, 

/c = / + 1, ...n. 



(2.9) 



In view of the mean value theorem 14 



f{ur{x),u^^\x), ^\x),x) - /(ur_i(a;),M^i\(x), ,u^J'_i^\x),x) 



n-l 



J2^ui'\x) f^(s){Ur-l{x),U^^li{x), ,uj."/\x),x) + 



=0 
n-l 



^ E ^Mr'Ha;) (^Mr^Ha;) (wr-i (a;), Mj.i\(x), ,mJ."/^(x),x), 



s,t=0 



where ) lies between ul^\x) and ^^.^(x). Now Eq. p.7| can be written as 



n-l 



L(")5m,+i(x) - E-^^lla;) A(«)K(x),u«(^), ,u("-i)(x),x) 



n-l 



^ E '^«r'n^)'^Wr*Ha;)/„(»)„w(Mr-i(a;),4-i(a;), , u^J'_i^\x) , x) . 



(2.10) 



(2.11) 



Denoting as the Greens function, which is the inverse of the following differential 



operator and incorporates linearized boundary conditions 2.5 and 2.6 



Lin) = Lin) .Y.U)Mx)M^\x), Mr'\x).x) — , (2.12) 

one can express the solution for the difference function 5ur+i as 

5Ur+l{x) = 

\ f' Gl^'\x,y)J2 H'\yM\y) U^ui^^i^^^^^ ,ut-i'\y),y)dy. (2.13) 

The functions 5u!;f\y)Su!ff\y) could be taken outside of the sign of the integral at some point 
y = X belonging to the interval, so one obtains 

-| 1 

5ur+iix) = -J2 5ui'\x)5ui'\x)Mstix). (2.14) 

2 s,t=0 

where Msr{x) equals 

M,,(x) = rG(")(x,y)/„(.)„w(^,_i(2/),wil^i(y), ,ut-i'\y),y)dy (2.15) 

Jo 

If Mst{x) is a strictly positive (negative) matrix for all x in the interval, then Sur+i{x) 
will be positive (negative), and the monotonic convergence from below (above) results. 
Obviously, from Eq. |2.13| follows 



\5Ur+lix)\ < kr{x)\\5Ur\f (2.16) 

where kr is given by 

kr{x) = - |G(")(x,y)| Y: \U)ui^){Ur-l{y),U^r\iy), , U^"/^ (?/) , ?/) M^/ (2.17) 

'^•^0 s,t=0 

and ||(5ur.|| is a maximal value of any of on the interval (0,b). 

Since Eq. |2.16| is correct for any x on the interval (0,b), it is correct also for some x = x 



where |5ur+i(x)| reaches its maximum value One therefore has 

\\SUr+l\\\ < krix)\\SUr\f (2.18) 



Assuming the boundness of the integrand in expression |2.17| , that is the existence of the 
bounding function F{x) such that integrand at a; = a; and at any y is less or equal to F{y), 
one finally has 

\\5Ur+l\\\ < k\\5ur\f , (2.19) 

where 

k = / F{x)dx . (2.20) 







The linearized boundary conditions |2.5| and p^.6| are obtained from exact boundary con 



ditions |2.2| and |2.3| by using the mean value theorem Eq. |2.1(J| and neglecting the quadratic 



terms, so that the error in using linearized boundary conditions vis-a-vis the exact ones 
is quadratic in the difference between the exact and linearized solutions. The maximum 
difference between boundary conditions |2]^ and ^]6] corresponding to two subsequent quasi- 
linear iterations is therefore quadratic in In view of this result and of Eq. |2.19| , the 



difference between the subsequent iterative solutions of Eq.|2.4| with boundary conditions 2.5 



and ^]6| decreases quadratically with each iteration. In a similar way, one can show [|l| that 
the difference Aur+i{x) = u{x) — Ur{x) between the exact solution and the r-th iteration is 
decreasing quadratically as well: 

||Am,.+i||| < k\\Aur.\f . (2.21) 

A simple induction of Eq. ^.19| shows |^ that 6un+i{x) for an arbitrary / < r satisfies 
the inequahty 

\\5ur+i\\<ik\\5ui+i\\f''-'/k, (2.22) 
or, for / = 0, we can relate the + 1 th order result to the 1st iterate by 

\\Sun+i\\ < ((^-ll^^i, 11)^7^- (2.23) 

The convergence depends therefore on the quantity qi = k\\ui — Uo\\, where, as we have 
mentioned earlier, the zeroth iteration Uo{x) is chosen from physical and mathematical 
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considerations. Usually it is advantageous (see discussion below) that uo{x) would satisfy 
at least one of the boundary conditions. From Eq. ( |2.22| ) it follows, however, that for 



convergence it is sufficient that just one of the quantities qm = k\\6um\\ is small enough. 
Consequently, one can always hope that even if the first convergent coefficient qi is 
large, a well chosen initial approximation uq results in the smallness of at least one of the 
convergence coefficients qm, ttl > 1, which then enables a rapid convergence of the iteration 
series for r > m. It is important to stress that in view of the quadratic convergence of the 
QLM method, the difference HA-Ur+iH between the exact solution and the QLM iteration 
always converges to zero if the difference 6ur+i{x) between two subsequent QLM iterations 
becomes infinitesimally small. 

Indeed, if 6ur{x) is close to zero, it means, since 6ur+i{x) = Aur{x) — Aur+i{x) that 
Aur{x) = Aur+i{x) OT Qr = Qr+1 whcrc Qr = k\\Aur\\. When one assumes the possi- 
bility that Qr and Qr+i could be not small, one could conclude that the iteration process 
"stagnates" , which means convergence to the wrong answer or no convergence at all. 

However, such a conclusion is wrong since Eq. p^.21| , which can be written as Qr+i < 



Ql, for Qr < I (this last inequality, starting from some r is a necessary condition of the 
convergence) could be not satisfied unless both ||(5r+i|| and ||Qr|| equal to zero. This proves 
that stagnation of the iteration process is impossible and convergence of H^Mf+iH to zero 
automatically leads to convergence of the QLM iteration sequence to the exact solution. 
Hence the QLM assures not only convergence, but also convergence to the correct solution. 

Another corollary of this iteration process is that if the solution and its derivatives are 
continuous functions of x, the convergence of the QLM in the whole region will follow. 
Indeed, even if the zero iteration uo{x) is chosen not to satisfy the boundary conditions, 
the next iteration Ui{x), being a solution of a linear equation with linearized boundary 



conditions p.5| and |2.6| , will automatically satisfy the exact boundary conditions |2]^ and | 
at least up to the second order in difference 6ui at the boundaries. This means that the 
difference between the exact and first QLM iterations at some intervals near the boundaries 
will be small, so that the QLM iterations in this interval would converge. Because the 
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subsequent values of k Sumix), m > 2 became much smaller for this interval, in view of 
assumed continuity of the solution and its derivatives these differences will also be small at 
the neighboring intervals. The subsequent iterations will extend the convergence to the next 
neighboring intervals and so on, until the convergence in the whole region will be reached. 
The predicted trend is therefore that the QLM yields rapid convergence starting at the 
regions where the boundary conditions are imposed and then spreading from there to all 
other regions. 

An additional important corollary is that, in view of Eq. |2.22| , once the quasilinear 



iteration sequence starts to converge, it will continue to do so, unlike the perturbation 
expansion, which is often given by an asymptotic series and therefore converges only up to 
a certain order and diverges thereafter. 

Based on this summary of the QLM, one can deduce the following important features of 
the quasilinearization method: 

i) The method approximates the solution of nonlinear differential equations by treating 
the nonlinear terms as a perturbation about the linear ones, and is not based, unlike 
perturbation theories, on the existence of some kind of small parameter. 

ii) The iterations converge uniformly and quadratically to the exact solution. In case of 
matrix M^t in Eq. |2.15| being strictly positive (negative) for all x in the interval, the 
convergence is also monotonic from below (above). 

iii) For rapid convergence it suffices that an initial guess for the zeroth iteration is suffi- 
ciently good to ensure the smallness of just one of the quantities g,. = k\u^.^i — u^H. If 
the solution and its derivatives are continuous, convergence follows from the fact that 
starting from the first iteration, all QLM iterations automatically satisfy the quasilin- 
earized boundary conditions |2.5| and p.6| . The convergence is extremely fast: if, for 
example, q\ is of the order of |, only 4 iterations are necessary to reach the accuracy 
of 8 digits, since (^)^" is of the order of (j^)^" ^ 
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iv) Convergence of H^Mr+iH to zero automatically leads to convergence of the QLM itera- 
tion sequence to the exact solution. 

v) Once the quasilinear iteration sequence at some interval starts to converge, it will 
always continue to do so. Unlike an asymptotic perturbation series, the quasilin- 
earization method yield the required precision once a successful initial guess generates 
convergence after a few steps. 



III. ANALYTICALLY SOLVABLE EXAMPLE: COMPARISON OF 
QUASILINEARIZATION APPROACH WITH EXACT SOLUTION AND WITH 

PERTURBATION THEORY. 

In order to investigate the applicability of the quasilinearization method and its conver- 
gence and numerical stability, let us start from a simple example of an analytically solvable 
nonlinear ordinary differential equation suggested in ref. p2| : 



u\r) = -g n"(r), u{0) = 1, (3.1) 

where the boundary condition at r = is also given and ' means differentiation in variable 
r. The exact solution to this problem is 

u{r) = {l + {n-l)gry^ (3.2) 

Since 

(l + x)^ = E ,r/V^ (3.3) 
g m!i [q + 1 — m) 



the expansion of the solution 3.2 in powers of g is given by 



^T(^)(g(n-1))'^ 

u{r) = E ^p/n 2 Y- ^"^ (3.4) 

n — m) 

U \n.— 1 / 

The convergence radius of the series |3.4] is i? = l/{g{n — 1)), which is inversely proportional 
to the extent n — 1 of the nonlinearity and to the value g of the perturbation parameter. 
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Now consider the quasilinearization approach to this equation, taking, for example, g = 1 
and n = Q. Here we consider Eq. |3.1| with a rather strong degree of nonhnearity. In this 
case, one can expect the convergence of the perturbation expansion only up to r < i? = |. 

The QLM procedure in the case where the nonlinear term depends only on the solution 
itself and not on its derivatives reduces to setting u'^_,_i(r) = f{uk) + {uk+i{r) —Uk{r)) fuiuk). 
Here / = —g u^{r) while its functional derivative fu equals to —g nu"'^^{r). The quasilin- 



earized equation p.4| for the {k + l)-th iteration for this case has therefore the following 
form: 

u[^^{r) + ngul~^{r)uk+i{r) = {n-l)gu'^{r), Uk+i{0) = I , (3.5) 

where Uk{r) is a previous iteration which is considered to be a known function. Let us choose 
as a zero iteration UQ{r) = 1 which satisfies the boundary condition Uq{0) = 1. 

The results of our QLM calculations with Eq. ^]5| are presented in Fig. |l] which displays 
the exact solution for the case of n = Q and g = 1, together with the first four QLM 
iterations. Convergence to the exact solution in Fig. |1| is monotonic from above as it should 
be as discussed in Section II and in Refs. [^0,^ due to fact that the second functional 
derivative —n{n — of the left-hand side of Eq. |3.1| for even n is strictly negative. 

The convergence starts at the boundary, exactly as expected from the discussion in section 
II, and expands with each iteration to larger values of the variable r. The difference between 
the exact solution and the sixth QLM iteration for all r in the range between zero and five 
where our calculations were performed is less than 10~^. Note that the QLM yields a solution 
beyond the convergence radius limit on the series solution of 1/5. 

IV. LANE-EMDEN EQUATION 

The Lane-Emden equation 

y"{r) + - y'{r) + y"(r) = 0, i/(0) = 1, y'{0) = (4.1) 
r 
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describes a variety of phenomena in theoretical physics and astrophysics, including aspects 
of stellar structure, the thermal history of a spherical cloud of gas, isothermal gas spheres, 
and thermionic currents, see ref. [17] and references therein. The parameter n defines an 
equation of state with its physically interesting range being < n < 5. The equation also 
appears in other contexts, e.g., in case of radiatively cooling, self gravitating gas clouds, in 
the mean-field treatment of a phase transition in critical absorption or in the modeling of 
clusters of galaxies. For n = 0, 1 and 5 the equation can be solved analytically. Setting 
y = ^ transforms the equation to a more convenient form without a first derivative: 

«"(^) + = 0' ^(0) = 0' = 1- (4-2) 

Let us consider this nonlinear equation for the physically interesting and analytically non- 



solvable case of n = 4. The quasilinearized form of equation 4.2 is 



. , , , urHr) _ n-1 



+ = ^ <{r), «fc+i(0) = 0, u'.^M = 1. (4.3) 

The simplest initial guess, satisfying the boundary conditions will be UQ{r) = r. Comparison 
of the quasilinear solutions corresponding to the first five iterations with the numerically 
computed exact solution are given in Fig. 0. The figure shows that the convergence to 
the exact solution is very fast. It starts, as in the example of the previous section, at the 
left boundary and spreads with each iteration to larger values of r as expected from the 
discussion in section 11. The difference between the exact solution and the eighth QLM 
iteration for all r in the range between zero and ten, where our calculations were performed, 
is less than 10^^^ . 

V. THOMAS-FERMI EQUATION 



The Thomas- Fermi equation p3| , pl | 



X u"{x) = u^{x), m(0) = 1, m(oo) = 0, (5.1) 



13 



is an equation for the electron density around the nucleus of the atom. The left hand side of 

the above equation equals zero for n < 0. The Thomas- Fermi equation is also very useful for 

calculating form-factors and for obtaining effective potentials which can be used as initial 

trial potentials in self-consistent field calculations.. It is also applicable to the study of 

nucleons in nuclei and electrons in metal. It is long known (see and references therein) 

that solution of this equation is very sensitive to a value of the first derivative at zero which 

insures smooth and monotonic decay from u{0) = 1 to u{oo) = as demanded by boundary 

conditions. By contrast, the computation is much simpler for the quasilinearized version 

of this equation. The QLM procedure in this case reduces to setting ^'^^.^(r) = f{uk) + 

— Uk{r))fu{uk), where / = " ^''^ and the functional derivative is fu = (3/2)^^, so 

that the QLM equation has a form: 

3 1 1 3 

Vx u'l^lix) - -Ul{x) Uk+i{x) = --Ul{x), Mfc+i(0) = 1, Mfc+i(oo) = 0, (5.2) 

which is easily solved by specifying directly the boundary condition at infinity without 
searching first for the proper value of the first derivative. The initial guess, satisfying the 
boundary condition at zero was chosen to be uo{x) = 1. The results of QLM calculations 
with Eq. are presented in Fig. ^ which displays the exact solution together with the 
first four QLM iterations. The convergence starts at the boundaries, exactly as expected 
from the discussion in section II, and expands with each iteration to a wider range of values 
of the variable x. The difference between the exact solution and the eighth QLM iteration 
for all X in the range between zero and forty where our calculations were performed is less 
than 10~'^. 

VI. CLASSICAL ANHARMONIC OSCILLATOR 

The classical anharmonic oscillator satisfies the nonlinear second-order Duffing equation 

u{t) +u{t) + g u^{t) =0 . (6.1) 
The typical initial conditions are 
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m(0) = 1, u{0) = 0. (6.2) 
The solution oscillates strongly and thus is more difficult to approximate. It is, for example, 



well known p2[ that the usual perturbative solution is valid only for times t small compared 



with i, so that for larger g the perturbative solution is adequate only on a small time 
interval. In contrast, the quasilinearization approach gives solution in the whole region also 
for large (7- values. 

The quasilinearized equation is 

Uk+i{t) + {l + 3 g ul{t))uk+i{t)-2 g ul{t) =0, Uk+i{0) = 1, Mfc+i(0) = 0. (6.3) 

The results of QLM calculations with Eq. for g = 3 are presented in Figs. ^ and |^. 
Fig. ^ displays the exact solution together with the QLM solutions for the ffist, second and 
fourth iterations while Fig. ^ shows comparison of exact solution with sixth, seventh and 
eighth QLM iterations. Again, the convergence starts at the left boundary as expected from 
the discussion in section II, and expands with each iteration to larger values of the variable 
t. The difference between the exact solution and the eleventh QLM iteration for all t in the 
range between zero and seven where our calculations were performed is less than 10"^'' . 



VII. BLASIUS EQUATION 

A third order nonlinear Blasius equation 

u"'{x) + u"{x)u{x) = 0, u(0) = u'{0) = , ?i'(oo) = 1 (7.1) 

describes the velocity proffie of the fluid in the boundary layer. The QLM procedure in this 
case is given by m'^Vi(^) = fi'^k, u'l) + {uk+i - Uk)fu{uk, u'l) + {u'l^^ - u'l)fu"{uk, u'l), where 
f{u,u') = —u"u, fu{u,u') = —u" and fu"{u,u') = —u. The quasilinearized version of the 
Blasius equation thus has a form 

Mfc+i(a;) + Uk{x)ul^^{x) + Uk+i{x)ul{x) - Uk{x)ul{x) = 0, 

Mfc+i(0) = 4+1(0) = ,4+1(00) = 1. (7.2) 
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The initial guess, satisfying the boundary condition for the derivative at zero was chosen 
to be Uq{x) = 1. The results of QLM calculations with Eq. [7.2| are presented in Fig. ^ which 
displays the exact solution together with the first QLM iteration. The convergence starts 
at the left boundary as follows from the discussion in section II, and expands with each 
iteration to larger values of the variable x. The difference between the exact solution and 
the fifth QLM iteration for all x in the range between zero and ten where our calculations 
were performed is less than 10~^^ . 

VIII. CONCLUSION 

Summing up, we formulated here the conditions under which the quadratic, uniform and 
often monotonic convergence of the quasilinearization method are valid. 

We have followed here the derivation outlined in ref. [|I|, which is not based, unlike the 
derivations in refs. on a smallness of the interval and on the boundness of the nonlinear 
term and its functional derivatives, the conditions which usually are not fulfilled in physical 
applications. 

In order to analyze and highlight the power and features of the quasilinearization method 
(QLM), in this work we have also made numerical computations on different ordinary second 
and third order nonlinear differential equations, applied in physics, such as the Blasius, 
Duffing, Lane-Emden and Thomas- Fermi equations and have compared the results obtained 
by the quasilinearization method with the exact solutions. Although all our examples deal 
only with linear boundary conditions, the nonlinear boundary conditions can be handled 
readily after their quasilinearization as explained in Section II. 

In all the examples considered in the paper the simplest initial guess was enough to 
produce a quadratic convergence. However, as it is well known for the Newton method on 
which quasilinearization method is based that the divergence could occur if the initial guess 
is too bad. In this case the modification of QLM based on the damped Newton method with 
relaxation factor /cite CB,RR may help. The price of such modification would be, like in 
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the damped Newton method only hnear convergence instead of the quadratic one. 
Our conclusions are as follows: 

The QLM treats nonlinear terms by a series of nonperturbative iterations and is not based 
on the existence of some kind of small parameter. At every iterative stage, the differential 
operator changes significantly to account for the nonlinearity, which is the major way that 
the QLM differs from other approximative techniques. As a result, as we see in all our 
examples, the QLM is able to handle large values of the coupling constant and any degree 
of the nonlinearity, unlike perturbation theory. Thus the QLM provides extremely accurate 
and numerically stable answers for a wide range of nonlinear physics problems. The QLM 
is also very easy to apply. 

In view of all this, since most equations of physics, from classical mechanics to quantum 
field theory, are either not linear or could be transformed into a nonlinear form, the quasi- 
linearization method appears to be extremely useful and in many cases more advantageous 
than the perturbation theory or its different modifications, like expansion in inverse powers 
of the coupling constant, the expansion, etc. 
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FIGURE CAPTIONS 



FIG. 1. Convergence of QLM iterations for the analytic example of section III and 
comparison with the exact solution. Thin solid, dot-dashed, short-dashed and dotted curves 
correspond to the first, second, third and fourth QLM iteration respectively, while the thick 
solid curve displays the exact solution. The convergence is monotonic from above as it 
should be according to the discussion in the text. The difference between the exact solution 
and the sixth QLM iteration for all r in the figure is less than 10~^. 

FIG. 2. Convergence of QLM iterations for the Lane-Emden equation and comparison 
with the numerically obtained exact solution. Thin sohd, dot-dashed, short-dashed, long- 
dashed and dotted curves correspond to the first, second, third, fourth and fifth QLM 
iteration, respectively, while the thick solid curve displays the exact solution. The difference 
between the exact solution and the eighth QLM iteration for all r in the figure is less than 
10-11 

FIG. 3. Convergence of QLM iterations for the Thomas-Fermi equation and comparison 
with the numerically obtained exact solution. Thin solid, dot-dashed, short-dashed and 
dotted curves correspond to the first, second, third and fourth QLM iteration, respectively, 
while the thick sohd curve displays the exact solution. The difference between the exact 
solution and the eighth QLM iteration for all x in the figure is less than 10"^ . 

FIG. 4. Convergence of the first few QLM iterations for the Duffing equation and compar- 
ison with the numerically obtained exact solution. The dotted curves on three consecutive 
graphs correspond to the first, second and fourth QLM iteration respectively, while the solid 
curve displays the exact solution. 

FIG. 5. Convergence of the higher QLM iterations for the Duffing equation and compari- 
son with the numerically obtained exact solution. The dotted curves on the three consecutive 
graphs correspond to the sixth, seventh and eighth QLM iteration respectively, while the 
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solid curve displays the exact solution. The difference between the exact solution and the 
eighth QLM iteration for all t in the figure is less than 10~^° . 

FIG. 6. Comparison of the first QLM iteration for the Blasius equation with the numer- 
ically obtained exact solution. The difi^erence between the exact solution and the fifth QLM 
iteration for all x in the figure is less than 10~^° . 
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FIG. 1. Convergence of QLM iterations for the analytic example of section III and comparison 

with the exact solution. Thin solid, dot-dashed, short-dashed and dotted curves correspond to 

the first, second, third and fourth QLM iteration respectively, while the thick solid curve displays 

the exact solution. The convergence is monotonic from above as it should be according to the 

discussion in the text. The difference between the exact solution and the sixth QLM iteration for 

all r in the figure is less than 10~^. 



u 




r 



2 4 6 8 10 



FIG. 2. Convergence of QLM iterations for the Lane-Emden equation and comparison with the 
numerically obtained exact solution. Thin solid, dot-dashed, short-dashed, long-dashed and dotted 
curves correspond to the first, second, third, fourth and fifth QLM iteration, respectively, while 
the thick solid curve displays the exact solution. The difference between the exact solution and 
the eighth QLM iteration for all r in the figure is less than 10~^^ . 
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FIG. 3. Convergence of QLM iterations for the Thomas-Fermi equation and comparison with 
the numerically obtained exact solution. Thin solid, dot-dashed, short-dashed and dotted curves 
correspond to the first, second, third and fourth QLM iteration, respectively, while the thick solid 
curve displays the exact solution. The difference between the exact solution and the eighth QLM 
iteration for all x in the figure is less than 10~^ . 
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FIG. 4. Convergence of the first few QLM iterations for the Duffing equation and comparison 
with the numerically obtained exact solution. The dotted curves on three consecutive graphs 
correspond to the first, second and fourth QLM iteration respectively, while the solid curve displays 
the exact solution. 
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FIG. 5. Convergence of the higher QLM iterations for the Duffing equation and comparison 
with the numerically obtained exact solution. The dotted curves on the three consecutive graphs 
correspond to the sixth, seventh and eighth QLM iteration respectively, while the solid curve 
displays the exact solution. The difference between the exact solution and the eighth QLM iteration 
for all t in the figure is less than lO"-*^^ . 
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FIG. 6. Comparison of the first QLM iteration for the Blasius equation with the numerically 

obtained exact solution. The difference between the exact solution and the fifth QLM iteration for 
all X in the figure is less than 10~^° . 
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